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Abstract 

Researchers  in  many  fields  often  need  to  quantify  the  similarity  between  images 
using  metrics  that  measure  qualities  of  interest  in  a  robust  quantitative  manner.  We 
present  here  the  concept  of  image  dimension  reduction  through  characteristic  shape  se¬ 
quences.  We  formulate  the  problem  as  a  nonlinear  optimization  program  and  demon¬ 
strate  the  solution  on  a  test  problem  of  extracting  maximal  area  ellipses  from  two- 
dimensional  image  data.  To  solve  the  problem  numerically,  we  augment  the  class  of 
mesh  adaptive  direct  search  (MADS)  algorithms  with  a  filter,  so  as  to  allow  infeasible 
starting  points  and  to  achieve  better  local  solutions.  Results  here  show  that  the  MADS 
filter  algorithm  is  successful  in  the  test  problem  of  finding  good  characteristic  ellipse 
solutions  from  simple  but  noisy  images. 

Key  words:  Optimization,  mesh  adaptive  direct  search  algorithms  (MADS), 
hlter  methods,  nonsmooth  optimization,  image  metrics 


1  Introduction 

In  this  paper  we  present  a  method  for  property  preserving  dimension  reduction  in  images 
that  can  be  used  to  construct  metrics  for  quantifying  similarity  between  two  images  in  a 
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high-dimensional  space.  This  should  not  be  confused  with  the  related  concept  of  determin¬ 
ing  image  measure  properties,  which  is  a  function  of  only  one  image.  Examples  of  image 
metrics  include  Lp  distances,  correlation  functions,  and  image  warping  distances.  The  latter 
capitalizes  on  image  registration  techniques  common  in  medical  imaging  applications  [20]. 
Image  metrics  are  also  essential  for  data  retrieval  applications 

The  most  useful  and  versatile  image  comparison  methods  will  share  significant  proper¬ 
ties.  A  key  idea  is  that  a  metric  must  quantitatively  measure  qualities  of  interest.  Equally 
important  is  that  a  metric  must  ignore  image  differences  that  do  not  matter  in  the  opinion 
of  an  expert.  One  such  difference  is  noise,  a  common  and  understandable  theme.  However, 
robustness  in  the  presence  of  non-noise  features,  such  as  distortions,  occlusions,  and  missing 
data,  is  also  important.  Metrics  are  especially  useful  if  they  are  automated  and  rely  as  little 
as  possible  on  expert  intervention. 

Our  consideration  is  the  characterization  of  an  image  through  extraction  of  geometric 
properties  of  image  intensity  level  sets  and  a  subsequent  metric  comparison.  The  problem 
that  inspired  our  task  is  parameter  estimation  and  hydrodynamic  code  validation  studies 
using  fluid-flow  experimental  data.  The  essential  question  is:  How  well  does  a  given  choice  of 
hydrodynamic  parameters  and  code  implementation  predict  what  is  recorded  in  experimental 
images?  The  related  question  is:  Can  we  dehne  a  metric  in  the  space  of  images  that  can 
help  us  make  better  parameter  choices  and  measure  the  quality  of  the  simulation  code? 
Experimental  fluid  images  are  typically  characterized  by  an  evolving  object  or  region  of 
interest  which  may  have  a  simple  geometry  such  as  a  wavefront,  or  a  more  complex  character 
such  as  the  onset  of  turbulent  flow.  In  either  case  we  can  expect  to  capitalize  on  geometric 
property  extraction.  For  other  methods  that  focus  on  geometric  object  properties,  see  m 

In  this  paper,  we  consider  a  method  of  characteristic  shape  sequences  as  a  reduced  di¬ 
mension  representation  of  images.  The  idea  is  to  optimally  map  each  intensity  level  set  of  the 
image  onto  a  parameterized  family  of  shapes.  If  the  shape  family  is  characteristic  of  image 
features  of  interest,  then  the  shape  parameter  sequences  dehne  a  low-dimensional  represen¬ 
tation  of  the  original  image.  Image  similarity  can  then  be  measured  by  sequence  similarity 
metrics.  This  simple  idea  may  be  useful  for  certain  classes  of  huid  dynamics  experiments  and 
medical  imaging  applications  for  which  characteristic  shapes  can  be  of  primary  importance, 
background  details  are  relatively  simple,  and  level  sets  are  not  topologically  complicated. 

The  concept  of  metric  comparisons  on  reduced  dimension  measures  is  not  new.  The 
idea  is  to  use  an  appropriate  mapping  from  the  high-dimensional  image  space  to  a  low¬ 
dimensional  measure  space  that  preserves,  as  much  as  possible,  the  image  qualities  of  interest, 
so  that  metrics  can  then  be  applied  to  the  reduced  image  data.  Recent  applications  include 
speech  summarization  ra.  hyperspectral  analysis  [15],  and  experimental  design  [10].  There 
are  two  good  reasons  for  this  approach.  First,  it  is  possible,  in  principle,  to  formulate  a 
property-preserving  mapping  to  a  reduced  dimension  space  [9]  [mile].  Second,  we  can  take 
advantage  of  existing  sequence  similarity  algorithms  for  metric  comparison  (see  [21]  and 
included  references). 
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As  a  test  case,  we  consider  the  sequence  of  maximal  area  ellipses  interior  to  image  level 
sets.  Each  ellipse  is  dehned  by  hve  parameters  in  the  image  space,  and  each  optimal  el¬ 
lipse  is  computed  by  solving  a  constrained  nonlinear  optimization  problem.  Robustness  of 
the  solution  to  noise  and  experimental  uncertainty  is  enforced  by  allowing  the  ellipses  to 
contain  a  specihed  small  percentage  of  pixels  that  lie  outside  the  specihed  image  level  set. 
The  optimization  problem  is  solved  numerically  by  applying  a  mesh  adaptive  direct  search 
(MADS)  algorithm  |6]  with  a  hlter  |5].  The  hlter  allows  intermediate  solutions  that  violate 
constraints  in  order  to  provide  a  more  robust  global  search  of  the  parameter  space. 

The  paper  is  outlined  as  follows.  In  Section]^ we  describe  and  formulate  the  problem  in 
detail.  In  Section]^ we  introduce  a  hybrid  mesh  adaptive  direct  search  hlter  (hlter-MADS) 
algorithm.  Computational  results  obtained  by  the  algorithm  are  presented  and  analyzed  in 
Section  followed  by  some  concluding  remarks  in  Section 


2  Problem  Formulation 

In  formalizing  the  mathematical  description  of  the  problem,  we  begin  with  a  pixilated  rect¬ 
angular  image  11  =  {ztj  =  z{xi,yj)  :  i  =  1,2, M,j  =  1,2, . . . ,  N},  where  {{xi,yj)  : 
i  =  1,  2, ... ,  M,j  =  1,2, ... ,  N}  C  are  the  coordinates  of  the  pixels  and  2;  =  z{x,  y) 
is  a  real- valued  function  that  measures  intensity  at  location  {x,y).  An  example  is  given 
in  Figure  [T|  where  intensity  values  are  diherentiated  by  color.  Without  loss  of  generality. 


Figure  1:  Test  Image  for  the  maximal  area  ellipse  problem 

we  assume  that  the  image  lies  in  the  hrst  quadrant  of  a  rectangular  coordinate  system  and 
that  the  coordinate  {xi,yi)  coincides  with  the  origin.  Thus,  for  each  pixel  {xi,yj)  G  11, 
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i  =  1,2, . . . ,  M,  j  =  1,2, . . . ,  N,  the  associated  pixel  indices  (i,j)  can  be  determined  by 
dividing  Xi  and  i/j  by  the  pixel  width  and  height  ph,  respectively. 

An  ellipse  E  =  E(xc,  Pc,  cl,  b,  6)  C  can  be  represented  by  hve  variables:  the  location  of 
its  center  {xdic)  G  1^^,  the  lengths  of  its  axes  a  >  0  and  6  >  0  (measured  from  the  center), 
and  the  angle  of  its  rotation  6  G  [0,  |],  measured  counterclockwise  from  the  positive  x-axis. 
Its  area  is  given  by  A  =  nab. 

The  objective  is  to  hnd  the  maximum  area  ellipse,  subject  to  reasonable  bound  constraints 
on  each  variable,  for  which  a  specihed  percentage  p  G  [0, 100]  of  the  2;-values  interior  to  the 
ellipse  meet  a  specihed  intensity  threshold  Zq  and  for  which  a  specihed  p  G  [0, 100]  percent 
of  the  values  on  the  ellipse  boundary  meet  the  threshold  Zq.  Specihcally,  the  optimization 
problem  can  be  formulated  by 


subject  to 


A  = 

=  nab 

z(x 

,y) 

> 

Zq 

for  p%  of  1 

[x,y)  e 

n  n  E{xc, 

Pc,  CL,  b,  6) 

z(x 

,y) 

> 

Zd 
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Xc 

< 
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Ph 

< 

Pc 
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< 

a 

< 
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0 

< 

b 

< 
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0 

< 

e 

< 

TT 

2  ’ 

?/max  — 

Nph- 

Since  no 

data  exists  for  a 

point  lying 

(1) 


where  Xmax  =  Mp^  a 

image  boundary  (dehned  by  the  bound  constraints  in  ([^),  objective  function  values  there 
are  assumed  to  be  inhnite.  This  is  not  a  problem  for  any  of  the  algorithms  used  here, 
provided  that  the  initial  point  satishes  the  bound  constraints. 


3  Algorithm 

3.1  Mesh  Adaptive  Direct  Search  (MADS) 

The  problem  given  in  ([^  is  not  solvable  by  traditional  gradient-based  methods  because  no 
derivative  information  for  the  nonlinear  constraint  functions  is  available.  Thus  we  turn  to 
derivative-free  methods,  and  in  particular,  the  class  of  mesh  adaptive  direct  search  (MADS) 
algorithms.  MADS  is  an  extension  of  the  class  of  generalized  pattern  search  (GPS)  methods 
[niiisiES]  to  general  nonlinear  programming  (NLP)  problems  that  overcomes  limitations 
of  GPS  when  dealing  with  nonlinear  constraints  and  nonsmoothness. 

Lewis  and  Torczon  [19]  extended  GPS  to  nonlinear  constraints  by  applying  GPS  to  an 
augmented  Lagrangian  subproblem  based  on  the  method  proposed  in  [7|.  Their  approach 
ensures  convergence  for  twice  continuously  differentiable  objectives  and  constraints  to  a  KKT 
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point.  This  is  not  a  realistic  assumption  for  the  class  of  problems  considered  in  this  paper, 
and,  in  practice,  the  augmented  Lagrangian  approach  may  suffer  from  many  of  the  problems 
associated  with  penalty  function  methods,  such  as  the  effective  choices  of  penalty  parameters 
and  Lagrange  multipliers. 

Audet  and  Dennis  j5]  introduced  a  GPS  hlter  method  (hlter-GPS),  in  which  new  iterates 
are  generated  whenever  simple  decrease  in  the  objective  or  an  aggregate  constraint  violation 
function  is  achieved.  Filter  methods  were  hrst  introduced  by  Fletcher  and  Leyffer  [12]  as  a 
way  to  globalize  sequential  linear  programming  (SLP)  and  sequential  quadratic  programming 
(SQP)  without  using  a  penalty  function.  The  Audet-Dennis  GPS  hlter  approach  avoids  the 
pitfalls  of  penalty  parameters,  but  its  convergence  to  standard  hrst-order  stationary  points 
is  not  guaranteed  [5]  due  to  GPS  limitations. 

Because  of  the  weaknesses  inherent  to  hlter-GPS,  Audet  and  Dennis  [B]  more  recently 
introduced  the  class  of  MADS  algorithms.  Instead  of  limiting  local  exploration  to  a  hnite 
number  of  directions  (as  GPS  does),  MADS  systematically  generates  an  asymptotically  dense 
set  of  directions  in  the  limit.  Because  of  this  property,  convergence  to  both  hrst-order  [6] 
and  second-order  [2]  stationary  points  is  guaranteed  under  reasonable  assumptions  without 
the  need  for  a  hlter  but  with  the  requirement  that  the  starting  point  be  feasible. 

One  primary  advantage  of  MADS  is  that  it  can  handle  problems  with  general  “Yes/No” 
or  set  constraints.  However,  MADS  requires  the  initial  point  to  be  feasible,  while  hlter-GPS 
does  not.  In  this  paper,  we  describe  a  hlter-MADS  algorithm,  which  we  show  in  Section]^ 
to  be  more  advantageous  than  MADS  and  hlter-GPS  in  solving  ([^. 

The  target  class  of  problems  is  of  the  form, 

min  f{x) 

x^X 

subject  to  C{x)  <  0, 

where  X  C  M”,  /  :  M”  — >•  M  and  C  =  (Gi,  G2, . . . ,  Cm)  '■  K"'  — >■  IP’"-  When  using  hlter-GPS, 
we  must  assume  that  the  set  Y  is  a  linearly  constrained  region;  otherwise,  convergence  to 
stationary  points  is  not  guaranteed.  However,  no  such  restriction  is  needed  for  MADS.  The 
algorithms  described  here  are  actually  applied,  not  to  the  function  /,  but  to  the  barrier 
function  fx  =  f  +  V’v,  where  xjjx  is  an  indicator  function  dehned  as  zero  inside  X  and  -|-cx) 
otherwise.  Gonvergence  of  the  algorithm  is  based  on  the  smoothness  properties  of  /,  not  fx- 

MADS  is  a  class  of  algorithms  that  generates  a  sequence  of  iterates  with  nonincreasing 
objective  function  values.  Associated  with  this  sequence  is  a  set  of  no  directions,  D  C  M", 
that  positively  span  M”;  i.e.,  any  vector  in  M”  can  be  represented  as  a  nonnegative  linear 
combination  of  vectors  in  D.  Furthermore,  each  direction  dj  G  D,  j  =  1,  2, . . .  ,nE>,  must  be 
constructed  such  that  dj  =  Gzj,  where  G  G  is  a  hxed  nonsingular  generating  matrix 

and  Zj  G  is  an  integer  vector.  For  convenience,  the  set  D  is  also  viewed  as  a  real  n  x  nj:, 
matrix  whose  columns  are  the  elements  of  the  set. 

At  each  iteration  k,  points  are  evaluated  on  a  mesh  constructed  using  vectors  in  D  in  an 
attempt  to  hnd  a  point  with  a  function  value  lower  than  any  previously  evaluated  points. 
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which  we  term  an  improved  mesh  point.  The  current  mesh  is  defined  as 

Mk=  \J{x  +  A^Dz  :  z  e 

x&Sk 

where  Sk  is  the  set  of  points  where  the  objective  function  /  had  already  been  evaluated 
by  the  start  of  iteration  k  [Sq  is  the  set  of  initial  points),  and  is  the  current  mesh  size 
parameter. 

Each  iteration  of  MADS  has  two  main  steps:  an  optional  search  and  a  local  poll.  In 
the  SEARCH  step,  any  finite  number  of  mesh  points  (including  none)  are  evaluated.  The 
strategy  for  picking  such  points  is  left  entirely  to  the  user.  Common  choices  are  discussed  in 
|1]  or  |5],  for  example. 

If  the  SEARCH  fails  to  find  an  improved  mesh  point,  then  the  more  rigidly  defined  poll 
step  is  invoked.  In  this  step,  a  subset  of  mesh  points  near  the  incumbent  are  evaluated, 
where  proximity  is  based  on  the  poll  size  parameter  A^.  This  set  of  trial  points  is  called 
a  frame.  Less  general  than  the  frames  of  Coope  and  Price  [8],  the  MADS  frame  uses  the 
current  incumbent  solution  Xk  (the  frame  center)  and  the  poll  and  mesh  size  parameters, 

>  0  and  A™  >  0,  respectively,  and  a  positive  spanning  set  of  directions  Dk  C  D. 
The  MADS  frame  at  iteration  k  is  more  formally  defined  |6]  to  be  the  set: 

Pk  =  {xk  +  A^d  :  d  G  Dk}  C  Mk, 

where  Dk  is  a  positive  spanning  set  such  that  0  ^  Dk,  and  for  each  d  E  Dk, 

•  d  can  be  written  as  a  nonnegative  integer  combination  of  the  directions  in  D  :  d  =  Du 

for  some  vector  u  G  that  may  depend  on  the  iteration  number  fc; 

•  the  distance  from  the  frame  center  Xk  to  a  frame  point  Xk  +  Affd  G  Pk  is  bounded  by 
a  constant  times  the  poll  size  parameter:  A™||(i||  <  A^max{||d'||  :  d'  G  D}; 

•  limits  (as  defined  in  Coope  and  Price  [S])  of  the  normalized  sets  Dk  are  positive  span¬ 
ning  sets. 

If  either  the  search  or  poll  step  is  successful  in  finding  an  improved  mesh  point,  then 
that  point  becomes  the  incumbent,  and  the  mesh  is  retained  or  coarsened.  If  both  the 
SEARCH  and  POLL  steps  fail  to  find  an  improved  mesh  point,  then  the  current  iterate  is 
declared  a  mesh  local  optimizer,  the  frame  is  called  a  minimal  frame,  and  the  frame  center 
Xk  is  a  minimal  frame  center.  In  this  case,  the  current  iterate  Xk+i  =  Xk  is  retained,  and  the 
mesh  is  refined.  Coarsening  and  refining  of  the  mesh  are  performed  according  to  the  rule, 

r""=A™  (2) 


\m 

^fc+l 


where  Wk 


G 


{0, 1, . . . ,  w"*"}  if  an  improved  mesh  point  is  found 
{  — 1,  — 2, . . . ,  ta”}  otherwise. 


(3) 


November  22,  2006 


7 


where  r  >  1  is  rational,  and  w  <  —  1  and  w~^  >  0  are  integers.  The  poll  size  parameter  is 
also  updated  such  that  A™  <  for  all  k,  but  also  such  that 

lim  A™  =  0  if  and  only  if  lim  A^  =  0,  (4) 

k£K  k£K 


where  K  denotes  any  subsequence  of  mesh  local  optimizers.  Under  this  construction,  GPS 
is  the  specihc  implementation  of  MADS  for  which  A^  =  A|  for  all  k. 

Figures]^ and give  examples  of  a  GPS  and  MADS  frame,  respectively  [H].  As  shown 


Ar  =  A'  =  1 


Ar 


1 

4 


Figure  2:  Example  of  GPS  frames  Pk  =  {xk  +  A™(i  :  d  G  Dk}  =  for  different 

values  of  A™  =  A^.  In  all  three  hgures,  the  mesh  Mk  is  the  intersection  of  all  lines. 


Figure  3:  Example  of  MADS  frames  Pk  =  {xk  +  A™d  :  d  G  Dk}  =  {p^,p‘^,p^}  for  different 
values  of  A™  and  A|.  In  all  three  hgures,  the  mesh  Mk  is  the  intersection  of  all  lines. 


November  22,  2006 


8 


in  Figure]^  the  number  points  eligible  to  be  chosen  for  a  MADS  poll  set  grows  unbounded 
as  the  iteration  sequence  progresses.  In  fact,  the  total  of  all  poll  directions  can  be  chosen 
to  become  asymptotically  dense.  A  specific  implementation  of  MADS  was  devised  in  [H]  in 
which  the  sets  are  limited  in  size  but  chosen  randomly  at  each  iteration.  This  preserves 
the  denseness  of  poll  directions,  but  does  so  only  with  probability  one. 

3.2  Filters 

Filter-based  algorithms  attempt  to  minimize  both  the  objective  function  /  and  a  nonnegative 
aggregate  constraint  violation  function  h,  where  h  satishes  the  condition  that  h{x)  >  0  with 
h{x)  =  0  if  and  only  if  x  is  feasible.  Consistent  with  |S],  we  specihcally  choose 

h{x)  =  ||C'(a;)  +  ||2, 

where  C'(x)+  is  the  vector  of  constraint  violations  at  x;  i.e.,  for  i  =  1,2, ...  ,m,  Cj(x)+  = 
Ci{x)  if  Ci{x)  >  0;  otherwise,  Ci(x)+  =  0.  This  choice  of  h  inherits  whatever  smoothness 
properties  C  possesses  [5]. 

Also  consistent  with  [5],  we  dehne  a  second  constraint  violation  function,  hx  =  h  +  iIjx, 
where  the  indicator  function  tfjx  is  dehned  exactly  as  before.  Thus,  for  x  ^  X,  hx{x)  =  oo, 
and  the  function  h  is  not  evaluated.  Convergence  properties  with  respect  to  h  at  the  limit 
point  of  the  algorithm  depend  on  the  local  smoothness  of  h  and  not  of  hx  [3]. 

A  hlter  is  based  on  the  idea  of  dominance,  which  is  defined  with  respect  to  /  and  h 
below  [3].  A  formal  hlter  dehnition  follows. 

Definition  3.1  A  point  x  eMA  is  said  to  dominate  y  G  M"',  written  x  -<  y,  if  f{x)  <  f{y) 
and  hx{x)  <  hxiy)  with  either  f{x)  <  f{y)  or  hx{x)  <  hxiy)- 

Definition  3.2  A  hlter  X  is  a  finite  set  of  points  in  the  domain  of  f  and  h  such  that  no 
pair  of  points  x  and  y  in  the  set  have  the  relation  x  -<  y. 

Two  additional  restrictions  are  placed  on  the  hlter  JF.  First,  a  bound  h^ax  is  set  on 
aggregate  constraint  violation  so  that  each  point  y  E  tF  satishes  hxiy)  <  hmax-  Secondly,  to 
be  consistent  with  only  infeasible  points  are  included  in  the  hlter,  and  feasible  points 

are  tracked  separately.  Using  these  two  restrictions,  we  include  the  following  terminology 

0: 

Definition  3.3  A  point  x  is  said  to  be  hltered  by  a  filter  T  if  any  of  the  following  properties 
hold: 


1-  y  A  X  for  some  y  E  iF , 
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A  GENERAL  MADS  FILTER  ALGORITHM 

1.  Initialization:  Let  Jq  be  a  set  of  initial  evaluated  points  in  X.  Include  in  the 
filter  Ao  all  undominated  points  in  Jq.  Set  /imax  >  hx{xo)  for  some  xq  G  Aq.  Let 
Aq  >  0,  and  set  D,  G,  r,  w~  and  satisfy  the  requirements  given  above.  Set  the 
iteration  counter  <—  0. 

2.  Choose  frame  center  pk  G  C  X. 

3.  Search  and  poll  step:  Perform  the  search  and  possibly  the  poll  steps  (or 
only  part  of  them)  until  an  unfiltered  mesh  point  Xk+i  G  is  found. 

•  Optional  search:  Evaluate  fx  and  hx  on  a  finite  subset  of  trial  points  on 
the  mesh  M^. 

•  Local  poll:  Evaluate  fx  and  hx  on  the  frame  Pk  centered  at  pk- 

4.  Parameter  update:  Update  according  to  and  so  that  0  is 

satisfied.  Update  filter  Pk+i  foi'  tbe  new  frame  center. 

Set  A;  <—  fc  +  1  and  go  back  to  Step  2. 

Figure  4:  A  general  MADS  filter  algorithm 


2.  hx{xj)  P 

3.  hx{x)  =  0  and  f{x)  >  f^,  where  =  min{/(w)  :  w  E  P',  h{w)  =  0}. 

The  point  x  is  said  to  he  unfiltered  by  T  if  it  is  not  filtered  by  T . 

Note  that  one  consequence  of  this  dehnition  is  that,  \iy  E  T ,  then  any  other  point  x  E 
satisfying  fix)  =  /(y),  h{x)  =  h{y)  cannot  be  added  to  the  hlter. 

3.3  The  MADS  Filter  Algorithm 

In  adding  the  hlter  to  the  MADS  algorithm,  we  change  notation  slightly.  Instead  of  polling 
around  the  current  iterate  we  poll  around  the  frame  center  pk  E  {pf ,  p^},  where  pf  G  A  is 
the  best  feasible  point  found  so  far  (i.e.,  the  feasible  point  with  the  lowest  objective  function 
value),  and  p\,  E  X  is  the  least  infeasible  point  (i.e.,  the  infeasible  point  with  the  lowest  value 
of  h).  If  no  feasible  point  has  been  found,  then  pk  =  p{  is  chosen.  Note  that,  by  construction, 
both  pf  and  p^,  and  thus  p^,  will  always  satisfy  the  constraints  that  dehne  X.  A  successful 
iteration  is  one  that  hnds  an  unhltered  point  is  found,  in  which  case,  it  is  added  to  the  hlter. 
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The  filter- MADS  algorithm  is  described  in  detail  in  Figure]^  Its  convergence  properties 
are  not  presented  here,  but  follow  directly  from  those  of  the  hlter-GPS  [5]  and  MADS  [6] 
algorithms. 


4  Numerical  Implementation 

To  evaluate  the  effectiveness  of  the  hlter-MADS  algorithm  on  the  maximal  area  ellipse 
problem,  we  tested  it,  together  with  the  hlter-GPS  and  MADS  algorithms,  using  the  NO- 
MADm  [1]  software  package.  The  test  image  (with  random  noise  already  added)  is  the  one 
illustrated  in  Figure  [T]  It  consists  of  a  256  by  260  matrix,  where  each  entry  corresponds 
to  the  2;-value  of  a  pixel,  ranging  in  value  from  —0.0637  to  1.0573.  Each  pixel  has  a  height 
of  ph  =  0.0125  and  a  width  of  =  0.01168.  The  2;-values  are  interpreted  as  the  discrete 
pixel-center  sampling  of  the  unknown  continuous  image. 

For  the  hlter-GPS  and  hlter-MADS  algorithms,  the  domain  X  is  dehned  by  the  hve 
bound  constraints  in  ([^,  while  the  other  two  nonlinear  constraints  are  treated  by  the  hlter. 
For  MADS,  X  is  dehned  by  all  of  the  constraints.  Recall  that  constraints  that  dehne  X  are 
treated  by  the  barrier  approach,  in  which  points  lying  outside  of  X  are  discarded  without 
being  evaluated. 


4.1  Initial  Iterates 


We  developed  two  diherent  methods  of  selecting  the  initial  iterate.  The  algorithm  is  able  to 
use  either  one,  or  a  combination  of  the  two.  In  the  hrst  approach,  we  set  the  initial  iterate 
to  be  a  modihed  version  of  the  inscribed  ellipse  of  the  image.  This  ellipse,  which  usually 
violates  the  hrst  two  nonlinear  constraints  in  ([^,  is  given  by 


Eq  =  E{xc,yc,a,b,9) 


E  ( 


:yn 


-a^max  -  ‘^PW, 


(5) 


In  the  second  approach,  we  compute  local  maxima  of  the  signed  distance  function  asso¬ 
ciated  the  level  set  {{x,y)  :  z{x,y)  >  zq}.  The  signed  distance  function  is  computed  using 
the  fast  marching  method  |22].  Each  local  maximizer  then  becomes  the  center  of  a  circle 
whose  radius  extends  90%  of  the  distance  to  the  boundary  of  the  level  set.  Each  circle  is 
then  given  an  orientation  so  that  the  diherence  in  intensity  between  the  two  diameter  end¬ 
points  is  maximal.  A  list  of  these  potential  initial  iterates  is  compiled,  and  any  circle  that 
signihcantly  overlaps  with  another  is  removed  from  the  list.  The  algorithm  then  evaluates 
each  ellipse  from  the  list  as  an  initial  iterate. 

The  hrst  approach  often  starts  infeasible,  but  with  an  ellipse  of  large  area.  This  often 
leads  to  better  solutions,  but  achieving  feasibility  can  be  a  challenge,  especially  for  high  values 
of  Zq.  On  the  other  hand,  the  second  approach  is  more  robust  in  Ending  a  feasible  point. 
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since  it  begins  with  smaller  ellipses,  but  it  is  also  more  likely  to  converge  to  a  local  solution. 
By  combining  the  two  methods,  a  third  approach  is  available,  in  which  the  algorithm  begins 
with  the  initial  point  Eq  and  then  uses  the  second  approach  if  the  algorithm  fails  to  hnd  a 
feasible  iterate  from  Eq. 


4.2  Test  Results 


We  solved  the  problem  numerically  using  the  two  different  initial  points,  as  described  in 
Section  44  When  using  the  hlter,  polling  at  each  iteration  was  done  around  both  the  best 
feasible  point  and  the  least  infeasible  point.  Whenever  either  the  MADS  or  hlter-MADS 
algorithm  was  used,  10  replications  were  performed  to  handle  the  inherent  random  factor. 
Problem  parameter  values  were  set  to  zq  =  0.241,  p  =  95%,  and  p  =  75%.  Termination  of 
the  algorithm  was  set  to  occur  when  the  mesh  size  parameter  fell  below  10“®. 


For  both  tests,  we  compared  the  performance  of  three  algorithms;  namely,  hlter-GPS  (us¬ 
ing  the  standard  standard  2n  directions),  MADS,  and  hlter-MADS.  In  the  case  of  MADS, 
infeasible  trial  points  were  simply  discarded.  When  using  the  MADS  or  hlter-MADS  algo¬ 
rithm,  we  recorded  the  best  result  from  the  10  replications  along  with  appropriate  averages. 

Tables  and  provide  results  for  feasible  and  infeasible  starting  points,  respectively. 


Table  1:  Results  when  starting  from  a  feasible  point. 


Evaluations 

Area 

Xc 

Vc 

a 

b 

6 

hlter-GPS 

1208 

1.52018 

1.3059 

2.051 

0.5763 

0.8396 

0.9132 

MADS 

best 

average 

403 

426 

1.58541 

1.44494 

0.8581 

2.334 

0.689 

0.7325 

0.5741 

hlter-MADS 

best 

average 

3161 

1987 

1.72286 

1.54562 

0.9192 

2.2875 

0.6981 

0.7856 

0.5509 

When  starting  from  a  feasible  initial  point,  all  three  algorithms  converged  to  a  local  mini- 
mizer,  with  hlter-MADS  hnding  a  better  solution  than  MADS  or  hlter-GPS.  However,  when 
the  infeasible  initial  point  was  used,  there  was  a  signihcant  diherence  in  results  between 
the  three  methods.  Although  hlter-GPS  was  able  to  hnd  a  feasible  ellipse,  it  was  far  from 
optimal.  This  is  not  surprising  because  only  two  of  the  ten  standard  2n  directions  used  by 
hlter-GPS  (increasing  a  or  h)  improve  the  objective  function  value,  since  the  other  three 
variables,  Xc,  Uc,  and  9,  have  no  ahect  on  the  objective  function  value.  In  contrast,  hlter- 
MADS  showed  signihcant  improvement.  The  expanded  set  of  search  directions  in  MADS 
allowed  the  algorithm  to  hnd  a  more  eccentric  ellipse  (i.e.,  one  with  a  higher  eccentricity 
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Table  2:  Results  when  starting  from  an  infeasible  point,  Eq. 


Evaluations 

Area 

Xc 

Vc 

a 

b 

9 

hlter-GPS 

757 

1.24320 

1.0186 

2.3461 

.7098 

.5575 

0 

MADS 

best 

average 

81 

81 

failed 

failed 

1.5184 

1.6 

1.575 

1.485 

0 

hlter-MADS 

best 

average 

488 

561 

1.79746 

1.13941 

1.2059 

2.0375 

0.575 

0.995 

0.5 

e  =  1^1  —  / a?)  to  better  fit  the  constraints.  Without  a  hlter,  MADS  did  not  move  from 

the  infeasible  initial  point.  Since  MADS  sets  the  objective  function  value  of  all  infeasible 
points  to  — cx)  (for  a  maximization  problem),  it  moves  as  soon  as  it  hnds  a  feasible  point. 
However,  in  this  case,  it  did  not  hnd  one.  This  is  consistent  with  the  theory,  which  requires 
a  feasible  initial  point  to  ensure  convergence  to  a  stationary  point  [6]. 

In  Figure]^  three  images  displaying  the  level  set  of  all  z  >  Zq  are  shown.  Superimposed 
over  each  image  is  the  boundary  (and  axes)  of  the  best  ellipse  found  by  the  hlter-GPS, 
MADS,  and  hlter- MADS  algorithms,  respectively,  when  starting  from  the  infeasible  point 
Eq. 


Figure  5:  Resulting  ellipses  starting  from  the  infeasible  point,  Eq\  hlter-GPS,  MADS,  and 
hlter-MADS,  respectively 

Once  we  were  able  to  solve  this  problem  for  a  particular  value  of  Zq,  we  constructed 
shape  sequences  by  computing  optimal  ellipses  for  a  full  range  of  Zq  values.  In  particular, 
we  partitioned  the  range  of  z-values  into  100  evenly  spaced  nodes,  Zi,i  =  1,2, ... ,  100,  ran 
10  replications  of  the  hlter-MADS  algorithm  for  each  node  (i.e.,  setting  zq  =  Zi  in  ([^),  and 
recorded  the  best  and  average  optimal  ellipses  for  each  run.  For  these  runs,  we  we  used  as 
starting  points  both  Eq  and  the  initial  points  found  using  the  signed  distance  function. 
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Figure  shows  a  plot  the  optimal  values  of  x,  y,  a,  b  and  9  for  each  value  of  zq,  and 
Figure  shows  the  best  and  average  optimal  objective  function  values  versus  Zq.  The  se¬ 
quences  in  Figure  represent  the  characteristic  signatures  of  the  test  image  in  the  ellipse 
representation.  The  plots  do  not  include  the  few  cases  in  which  the  algorithm  failed  to  hnd 
a  feasible  ellipse.  In  addition,  other  vector  quantities,  such  as  eccentricity,  can  be  obtained 
from  the  basic  sequences. 

The  high-frequency  variation  in  the  ellipse  signatures  can  be  attributed  to  two  sources. 
First,  and  most  signihcant,  hlter-MADS  is  hnding  a  locally  optimal  solution  for  each  zo  on 
an  objective  function  surface  complicated  by  signihcant  noise.  Second,  in  the  presence  of 
noise,  the  parameters  will  show  a  signihcant  sensitivity  to  Zq,  even  if  the  global  solution  is 
attained  for  each  value  of  zo- 

The  sequences  of  Figure  |^can  be  compared  to  identically  obtained  sequences  on  other 
images  as  an  image  metric,  but  these  ideas  are  beyond  the  scope  of  this  paper. 


Figure  6:  Evolution  of  the  optimal  ellipse  as  values  of  Zq  increase 


5  Conclusions 

We  have  demonstrated  a  method  for  quantifying  reduced-dimension  measures  on  images 
which  can  be  used  for  further  exploration  as  image  comparison  metrics.  The  maximal  area 
ellipse  problem  explored  in  this  paper  contains  the  essential  features  of  this  class  of  problems: 
whole  image  characterization,  low-dimensional  characteristic  property  representation,  and 
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Figure  7:  Evolution  of  the  objective  function  as  values  of  zq  increase 

robustness  to  noise.  Other  potential  interesting  characteristic  shapes  include  models  of 
bones,  library  templates  for  image  classihcation,  simple  closed  curves  (rectangles,  ellipses, 
stars),  and  alphanumeric  characters. 

While  hlter-GPS  and  MADS  are  both  appropriate  methods  for  treating  the  general  class 
of  problems,  the  hybrid  hlter-MADS  approach  achieves  a  better  result,  primarily  due  to  the 
superior  convergence  properties  of  MADS  and  the  flexibility  of  hlter  methods  in  starting 
from  a  good  but  infeasible  point. 

The  hlter-MADS  solutions  do  show  signihcant  solution  variability.  This  shortcoming  can 
be  addressed,  as  it  was  in  this  study,  by  considering  the  best  of  several  runs.  Alternatively, 
some  image  preprocessing  to  reduce  noise  signature  is  a  reasonable  analysis  step. 
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